
** rr_results_for_replication: Runs the primary PA production/liability analysis, along with some data cleaning for the violations analysis.

clear all

** Set major file path

global MASTER "C:\Users\tpeis\Dropbox\Coal Bonding"

** Standard Header Stuff
capture log close
set more off

************************************************************************************
************************************************************************************

set scheme stmono2

cd "$MASTER\Data\Mine Level Production\PA\"

use "$MASTER\Data\Mine Level Production\PA\pa_osmre_cleaned_annual", clear


* Some key variable labels
label var b05 "Bituminous x Post-05"
label var u05 "Underground x Post-05"
label var b01 "Bituminous x Post-01"
label var treatment_05 "Post-05"
label var treatment_01 "Post-01"
label var bunderground "Bituminous x Underground"



* Sort and panel set
sort permid year
xtset permid year

* Generate mine age
sort permid year
by permid: gen age = _n
gen ln_age = ln(age)

* Generate minimum and maximum ages by permit
bysort permid: egen min_log_age = min(ln_age)
bysort permid: egen max_log_age = min(ln_age)


* Generate log liability (lct = log combined total)
gen lct = ln(combined_total)

* By mine instead of permit
bysort mine_id year: gen ct_mine = sum(combined_total)
gen lct_mine = ln(ct_mine)


* Generate "No Liability" variable equivalent to dropping out of liability data (merge_with_0 is the variable that determines the liability part)
sort permid year
gen no_liability = F.log_tons==. & F.merge_with_0!=2 & FF.log_tons==. & FF.merge_with_0!=2
bysort mine_id year: egen nol_mine = max(no_liability)


* Generate Table 1
gen tons_000 = NetTonsTonsaftermoisturea/1000
gen ct_000 = combined_total/1000

** Table 1
estpost tabstat ct_000 if year>=1996 & year<=2014, by(bituminous) statistics(mean sd count max)

** Table 2
* Log liability and log tons by mine type
ttest lct if year<=2001 & year>=1996 & surface==1, by(bituminous)
ttest log_tons if year<=2001 & year>=1996 & surface==1, by(bituminous)


ttest lct if year<=2001 & year>=1996 & surface==0, by(bituminous)
ttest log_tons if year<=2001 &  year>=1996 & surface==0, by(bituminous)




** Table 3 (liability for surface)
** Table 3 column 1
sort permid year
reghdfe lct b01 L.log_tons ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1, absorb(year permid) cluster(permid)

* Table 3 column 2
reghdfe lct b01 log_tons ln_age lppi ln_pop if year>=1996 & year<2014 & surface==1, absorb(year permid) cluster(permid)

* Table 3 Column 3
reghdfe lct b01 log_tons ln_age lppi ln_pop bituminous if year>=1996 & year<2014 & surface==1, absorb(year) cluster(permid)

* Table 3 column 4
bysort county year: egen county_mines = count(permid)
gen lncm = ln(county_mines)
* IV Regressions
sort permid year
ivreghdfe lct b01 ln_age lppi ln_pop (log_tons = L.lncm) if year>=1996 & year<2014 & surface==1, absorb(year permid) cluster(permid)

*ivreghdfe lct b01 ln_age lppi ln_pop (log_tons = L.lncm) if year>=1996 & year<2014 & surface==1, absorb(year company_id) cluster(permid)





** Table 4 (liability for underground)

* Table 4 column 1
reghdfe lct b01 L.log_tons ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==0, absorb(year permid) cluster(permid)

* Table 4 column 2
reghdfe lct b01 log_tons ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==0, absorb(year permid) cluster(permid)

* Table 4 column 3
reghdfe lct b01 ln_age lppi ln_pop log_tons bituminous if year>=1996 & year<2014 & surface==0, absorb(year) cluster(permid)

* Table 4 column 4 (IV)
bysort county year: egen county_u = sum(underground)
gen lncu = ln(county_u)

sort permid year
ivreghdfe lct b01 ln_age lppi ln_pop (log_tons = L.lncu) if year>=1996 & year<2014 & surface==0, absorb(year permid) cluster(permid)


** Table 5
* Generate dose treatment
bysort permid treatment_01: egen avg_liability_pre_post = mean(lct) if year>=1996

gen avg_liability_post = avg_liability_pre_post if treatment_01==1
bysort permid: egen avg_liability_post_full = max(avg_liability_post)

gen avg_liability_pre = avg_liability_pre_post if treatment_01==0
bysort permid: egen avg_liability_pre_full = max(avg_liability_pre)

gen dose = avg_liability_pre_full*bituminous*treatment_01


* Generate CEM Matching for ABS policy (doing both population and density for later robustness)
* Population
cem lct ln_pop if year==2000 & surface==1, treatment(bituminous) k2k

bysort permid: egen cem_matched_full = max(cem_matched) if surface==1

drop cem_matched

* Pop density
cem lct popdens if year==2000 & surface==1, treatment(bituminous) k2k

bysort permid: egen cem_matched_dens_surf = max(cem_matched) if surface==1

drop cem_matched


* Column 1
* Pop
reghdfe bond_miss_prod_0 dose bituminous ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_full==1, absorb(year) cluster(permid)


* Column 2
* Pop
reghdfe bond_miss_prod_0 dose ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_full==1, absorb(permid year) cluster(permid)
*reghdfe bond_miss_prod_0 b01 ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1, absorb(permid year) cluster(permid)

* Density



* Column 3
* Pop
reghdfe bond_miss_prod_0 dose ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_full==1, absorb(permid year company_id) cluster(permid)

* Column 4 - not matched
reghdfe bond_miss_prod_0 dose ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1, absorb(permid year) cluster(permid)


* Binary checks for ABS

* Column 5
* Pop
reghdfe bond_miss_prod_0 b01 ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_full==1, absorb(permid year) cluster(permid)

*Column 6
reghdfe bond_miss_prod_0 b01 ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1, absorb(permid year) cluster(permid)



** Table 6 in other file

** Table 7

* Column 1
reghdfe no_liability dose lppi ln_age ln_pop bituminous if year>=1996 & year<=2014 & surface==1, absorb(year) cluster(permid)

* Column 2
reghdfe no_liability dose lppi ln_age ln_pop if year>=1996 & year<=2014 & surface==1, absorb(year company_id) cluster(permid)

* Column 3
reghdfe no_liability dose lppi ln_age ln_pop bituminous if year>=1996 & year<=2014 & surface==1 & cem_matched_full==1, absorb(year) cluster(permid)

* Column 4
reghdfe no_liability b01 lppi ln_age ln_pop bituminous if year>=1996 & year<=2014 & surface==1, absorb(year) cluster(permid)

* Column 5
reghdfe no_liability b01 lppi ln_age ln_pop if year>=1996 & year<=2014 & surface==1, absorb(year company_id) cluster(permid)

* Column 6
reghdfe no_liability b01 lppi ln_age ln_pop bituminous if year>=1996 & year<=2014 & surface==1 & cem_matched_full==1, absorb(year) cluster(permid)


** Table 8 in other file

** Table 9 in other file

** Table 10 in other file

** Table 11

** CEM Matching for ABS - 1996

cem lct ln_pop if year==1996 & surface==1, treatment(bituminous) k2k

bysort permid: egen cem_matched_full_96 = max(cem_matched) if surface==1

drop cem_matched

*Row 1
reghdfe bond_miss_prod_0 dose ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_full==1, absorb(permid year) cluster(permid)

reghdfe bond_miss_prod_0 dose ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_full_96==1, absorb(permid year) cluster(permid)

* Row 2 in other file

** Table 12

* Column 1
reghdfe log_tons dose bituminous ln_age lppi ln_pop if surface==1 & year>=1996 & year<=2014 & cem_matched_full==1, absorb(year) cluster(permid)

* Column 2
reghdfe log_tons dose ln_age lppi ln_pop if surface==1 & year>=1996 & year<=2014 & cem_matched_full==1, absorb(year permid) cluster(permid)

* Column 3
reghdfe log_tons b01 bituminous ln_age lppi ln_pop if surface==1 & year>=1996 & year<=2014 & cem_matched_full==1, absorb(year) cluster(permid)

* Column 4
reghdfe log_tons b01 ln_age lppi ln_pop if surface==1 & year>=1996 & year<=2014 & cem_matched_full==1, absorb(year permid) cluster(permid)


** Table 13
* Try excluding companies that own both kinds of mines
bysort company_id: egen max_firm_surface = max(surface)
bysort company_id: egen max_firm_underground= max(underground)


* Row 1
sort permid year
reghdfe lct b01 L.log_tons ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & max_firm_underground==0, absorb(year permid) cluster(permid)


* Row 2
reghdfe  bond_miss_prod_0 dose ln_age lppi ln_pop if year>=1996 & year<=2014 &  surface==1 & max_firm_underground==0 & cem_matched_full==1, absorb(year permid) cluster(permid)


* Get combined variable
gen abandoned = NetTonsTonsaftermoisturea==0
replace abandoned = . if NetTonsTonsaftermoisturea==.
gen either_0 = abandoned==1 | bond_miss_prod_0==1

* Row 3
reghdfe either_0 dose ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_full==1 & max_firm_underground==0, absorb(permid year) cluster(permid)

* Row 4
reghdfe no_liability dose lppi ln_age ln_pop bituminous if year>=1996 & year<=2014 & surface==1 & max_firm_underground==0 & cem_matched_full==1, absorb(year) cluster(permid)

* Rows 5 and 6 in other file


** Table 14
* Try sone nnmatch robustness - just one match
sort permid year
gen llagtons = LL.log_tons
gen lllagtons = LLL.log_tons

* Idle
* Row 1
teffects nnmatch (bond_miss_prod_0 ln_age lppi ln_pop max_prod llagtons lllagtons year treatment_01) (b01) if surface==1 & year>=1996 & year<=2014, dmv vce(robust) nneighbor(1) atet tlevel(1) con(0)

* Row 2
teffects nnmatch (either_0 ln_age lppi ln_pop max_prod llagtons lllagtons year treatment_01) (b01) if surface==1 & year>=1996 & year<=2014, dmv vce(robust) nneighbor(1) atet tlevel(1) con(0)

* Row 3
teffects nnmatch (no_liability ln_age lppi max_prod ln_pop year treatment_01) (b01) if surface==1 & year>=1996 & year<=2014, dmv nneighbor(1) vce(robust) atet

* Row 4 in other file

** Table 15

* Row 1
reghdfe bond_miss_prod_0 dose ln_age lppi ln_pop if year>=1996 & year<=2010 & surface==1 & cem_matched_full==1, absorb(permid year) cluster(permid)

* Row 2
reghdfe no_liability dose lppi ln_age ln_pop bituminous if year>=1996 & year<=2010 & surface==1 & cem_matched_full==1, absorb(year company_id) cluster(permid)

* Row 3 in other file

** Table 16 - Density

* Column 1 Row 1
reghdfe bond_miss_prod_0 dose bituminous ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_dens_surf==1, absorb(year permid) cluster(permid)

* Column 2 Row 1
reghdfe bond_miss_prod_0 b01 ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_dens_surf==1, absorb(permid year) cluster(permid)

* Column 1 Row 2
reghdfe no_liability dose bituminous ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_dens_surf==1, absorb(year) cluster(permid)

* Column 2 Row 2
reghdfe no_liability b01 bituminous ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==1 & cem_matched_dens_surf==1, absorb(year) cluster(permid)

* Row 3 in other file


*** Appendix

** Table 17

* CEM matching for water claims
cem lct ln_pop if year==2004 & surface==0, treatment(bituminous) k2k

bysort permid: egen cem_matched_full_sub = max(cem_matched) if surface==0

drop cem_matched




* Column 1
reghdfe bond_miss_prod_0 b05 ln_age lppi ln_pop bituminous if year>=1996 & year<=2014 & surface==0, absorb(year) cluster(permid)

* Column 2
reghdfe bond_miss_prod_0 b05 ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==0, absorb(year permid) cluster(permid)

* Column 3
reghdfe bond_miss_prod_0 b05 ln_age lppi ln_pop if year>=1996 & year<=2018 & surface==0, absorb(year permid) cluster(permid)


* Matched
* Column 4
reghdfe bond_miss_prod_0 b05 ln_age lppi bituminous ln_pop if year>=1996 & year<=2014 & surface==0 & cem_matched_full_sub==1, absorb(year) cluster(permid)

* Column 5
reghdfe bond_miss_prod_0 b05 ln_age lppi ln_pop if year>=1996 & year<=2014 & surface==0 & cem_matched_full_sub==1, absorb(year permid) cluster(permid)

* Column 6
reghdfe bond_miss_prod_0 b05 ln_age lppi ln_pop if year>=1996 & year<=2018 & surface==0 & cem_matched_full_sub==1, absorb(year permid) cluster(permid)



** Table 18
* Column 1
reghdfe no_liability b05 ln_age lppi ln_pop bituminous if year>=1996 & surface==0 & year<=2014, absorb(year) cluster(permid)
* Column 2
reghdfe no_liability b05 ln_age lppi ln_pop if year>=1996 & surface==0 & year<=2014, absorb(year company_id) cluster(permid)

* Column 3
reghdfe no_liability b05 ln_age lppi ln_pop bituminous if year>=1996 & surface==0 & year<=2014, absorb(year) cluster(permid)

* Column 4
reghdfe no_liability b05 ln_age lppi ln_pop if year>=1996 & surface==0 & year<=2018, absorb(year company_id) cluster(permid)


** Tables 19 and 20 in other files



reghdfe no_liability b05 ln_age lppi bituminous if year>=1996 & surface==0 & year<=2018 & max_firm_surface==0, absorb(year) cluster(permid)
reghdfe no_liability b05 ln_age lppi bituminous if year>=1996 & surface==0 & year<=2018, absorb(year) cluster(permid)
reghdfe no_liability b05 ln_age lppi if year>=1996 & surface==0 & year<=2018, absorb(year company_id) cluster(permid)
reghdfe no_liability b05 ln_age lppi bituminous if year>=1996 & surface==0 & year<=2019, absorb(year) cluster(permid)


** Table 21


** Water claim with dose
bysort permid treatment_05: egen avg_liability_pre_post_water = mean(lct) if year>=1996

gen avg_liability_post_water = avg_liability_pre_post_water if treatment_05==1
bysort permid: egen avg_liability_post_water_full = max(avg_liability_post_water)

gen avg_liability_pre_water = avg_liability_pre_post_water if treatment_05==0
bysort permid: egen avg_liability_pre_water_full = max(avg_liability_pre_water)

gen dose_water = avg_liability_pre_water_full*bituminous*treatment_05


* Column 1
reghdfe log_tons dose_water ln_age lppi ln_pop if surface==0 & year>=1996 & year<=2014, absorb(year permid) cluster(permid)

* Column 2
reghdfe log_tons b05 bituminous ln_age lppi ln_pop if surface==0 & year>=1996 & year<=2014, absorb(year) cluster(permid)

* Column 3
reghdfe log_tons b05 ln_age lppi ln_pop if surface==0 & year>=1996 & year<=2014, absorb(year permid) cluster(permid)


** Table 22

* Column 1
reghdfe bond_miss_prod_0 dose_water ln_age lppi bituminous if year>=1996 & year<=2014 & surface==0, absorb(year) cluster(permid)

* Column 2
reghdfe bond_miss_prod_0 dose_water ln_age lppi if year>=1996 & year<=2014 & surface==0, absorb(year permid) cluster(permid)




*** Figures

** Figure 1

** PPI graph
bysort year bituminous: egen year_mean_ppi = mean(ppi)

twoway (connected year_mean_ppi year if year>=1996 & bituminous==1 & year<=2014) (connected year_mean_ppi year if year>=1996  & bituminous==0 & year<=2014), xtitle("Year") ytitle("Mean PPI") legend(lab(1 "Bituminous") lab(2 "Anthracite")) ///
saving(mean_ppi, replace)

graph export "$MASTER\Tables\pa_mean_ppi.pdf", replace




** Figure 2 available at https://www.eia.gov/coal/markets/includes/archive2.php

** Figure 3

gen ones = 1
bysort year: egen mean_prod = mean(tons_000)
bysort year: egen mean_N = sum(ones)

twoway (connected mean_prod year if year>1996 & year<2015), xtitle("Year") ytitle("Average Production (000 Tons)") ///
saving(prod_time, replace)
 
graph export "$MASTER\Tables\prod_time.pdf", replace


** Figure 4

* Generate relevant percentages
bysort year bituminous underground: egen yearly_0 = sum(abandoned)
bysort year bituminous underground: egen yearly_count = sum(ones)
bysort year bituminous underground: egen yearly_active_count = sum(ones*(1 - bond_miss_prod_0)*(1-abandoned))
bysort year bituminous underground: egen yearly_nol = sum(no_liability)
bysort year bituminous underground: egen yearly_count_full = sum(ones)

gen abandoned_pct = yearly_0/yearly_count

bysort year bituminous underground: egen yearly_idle = sum(bond_miss_prod_0)

gen idle_pct = yearly_idle/yearly_count

gen exit_pct = abandoned_pct + idle_pct

gen resolve_pct = yearly_nol/yearly_count_full

sort year

* Graphs
twoway (connected yearly_active_count year if year>=1996 & year<=2018 & surface==1 & bituminous==1 & cem_matched_full==1) (connected yearly_active_count year if year>=1996 & year<=2018 & underground==0 & bituminous==0 & cem_matched_full==1), xtitle("Year") ytitle("Total Active Surface Mines") legend(lab(1 "Bituminous") lab(2 "Anthracite")) ///
saving(prod_total, replace)
graph export "$MASTER\Tables\pa_active_surface.pdf", replace


twoway (connected yearly_active_count year if year>=1996 & underground==1 & bituminous==1) (connected yearly_active_count year if year>=1996 & underground==1 & bituminous==0), xtitle("Year") ytitle("Total Active Underground Mines") legend(lab(1 "Bituminous") lab(2 "Anthracite")) ///
saving(prod_total, replace)
graph export "$MASTER\Tables\pa_active_underground.pdf", replace




** Figure 5

twoway (connected exit_pct year if year>=1996 & year<=2014 & surface==1 & bituminous==1 & cem_matched_full==1) (connected exit_pct year if year>=1996 & year<=2014 & underground==0 & bituminous==0 & cem_matched_full==1), xtitle("Year") ytitle("Proportion Mines Not Operating") legend(lab(1 "Bituminous") lab(2 "Anthracite")) ///
saving(prod_total, replace)
graph export "$MASTER\Tables\pa_percent_0_surface.pdf", replace

twoway (connected exit_pct year if year>=1996 & year<=2018 & underground==1 & bituminous==1) (connected exit_pct year if year>=1996 & underground==1 & bituminous==0 & year<=2018), xtitle("Year") ytitle("Proportion Mines Not Operating") legend(lab(1 "Bituminous") lab(2 "Anthracite")) ///
saving(prod_total, replace)
graph export "$MASTER\Tables\pa_percent_0_underground.pdf", replace




** Figure 6
twoway (connected resolve_pct year if year>=1996 & year<=2018 & underground==0 & bituminous==1) (connected resolve_pct year if year>=1996 & underground==0 & bituminous==0 & year<=2018), xtitle("Year") ytitle("Mines Fully Resolved") legend(lab(1 "Bituminous") lab(2 "Anthracite")) ///
saving(prod_total, replace)
graph export "$MASTER\Tables\pa_nol_count_surface.pdf", replace


twoway (connected resolve_pct year if year>=1996 & underground==1 & bituminous==1 & year<=2018) (connected resolve_pct year if year>=1996 & underground==1 & bituminous==0 & year<=2018), xtitle("Year") ytitle("Mines Fully Resolved") legend(lab(1 "Bituminous") lab(2 "Anthracite")) ///
saving(prod_total, replace)
graph export "$MASTER\Tables\pa_nol_count_underground.pdf", replace


** Figure 7 in other file

** Figure 8


* Event study for ABS liability

gen ytreat_b = bituminous*(year - 2001)
gen pre_trend_abs = ytreat_b*(year<=2001)
replace ytreat_b = . if bituminous==0


eventdd lct log_tons ln_age ln_pop if year>=1996 & year<2014 & surface==1, timevar(ytreat_b) hdfe absorb(year permid) cluster(permid)
graph export "$MASTER\Tables\pa_abs_liability.pdf", replace
estat leads
* leads p-value of .882




* Figure 9
* Coefplot for "idle mines"
gen dose_es = avg_liability_pre_full*bituminous


reghdfe bond_miss_prod_0 c.dose_es##ib2000.year ln_age ln_pop if year>=1996 & year<=2014 & cem_matched_full==1, absorb(year permid) cluster(permid)

matrix list e(b)

*pretrends power 0.2, pre(1/4) post(5/18) omit
*pretrends, pre(1/4) post(5/18) omit slope(`r(slope)')

test 1996.year#c.dose_es 1997.year#c.dose_es 1998.year#c.dose_es 1999.year#c.dose_es 2000.year#c.dose_es
* p-value for all leads: .851
test 1997.year#c.dose_es 1998.year#c.dose_es 1999.year#c.dose_es
* p-value for middle leads: .951

label variable year ""
coefplot, keep(*c.dose_es) nolabels coeflabels(, truncate(4)) baselevels omitted vertical xline(5) yline(0) xlabel(, alternate angle(90) labsize(small)) xtitle("Year") ytitle("Coefficient Value")


graph export "$MASTER\Tables\pa_idle_01.pdf", replace


* Figure 10 in other file


** Appendix

* Figure 11

* Try distributed lag
gen ytreat = bituminous*underground*(year - 2005)
gen ytreat_pre = ytreat*(year<2005)
replace ytreat = . if bunderground==0




gen ytreat_s = surface*(year - 2001)

replace ytreat_s = . if surface==0

* Water claim idle event study
eventdd bond_miss_prod_0 ln_age ln_pop if year>=1996 & year<=2014 & surface==0, timevar(ytreat) hdfe absorb(permid year) cluster(permid) graph_op(saving("$MASTER\Tables\pa_bmp_ddd", replace))
graph export "$MASTER\Tables\pa_bmp_sub.pdf", replace
estat leads
* p value of .046



* Figure 12

* Water claim no liability event studies
eventdd no_liability ln_age ln_pop bituminous if year>=1996 & surface==0 & year<=2018, timevar(ytreat) hdfe absorb(year) cluster(permid) graph_op(saving("$MASTER\Tables\pa_nol_ddd", replace))
graph export "$MASTER\Tables\pa_nol_ddd.pdf", replace

eventdd no_liability ln_age ln_pop if year>=1996 & year<=2018 & surface==0 & max_firm_surface==0, timevar(ytreat) hdfe absorb(year) cluster(permid) graph_op(saving("$MASTER\Tables\pa_nol_ddd", replace))
graph export "$MASTER\Tables\pa_nol_ddd_nos.pdf", replace



** Figure 13 in other file



* Generate other file for violations analysis
cd "$MASTER\Data\Mine Level Production\PA\"

save "pa_for_violations.dta", replace

